slow switching in a population of delayed pulse-coupled oscillators 
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We show that peculiar collective dynamics called slow switching arises in a population of leaky 
C^*) ■ integrate-and-fire oscillators with delayed, all-to-all pulse-couplings. By considering the stability of 

cluster states and symmetry possessed by our model, we argue that saddle connections between a 
pair of the two-cluster states are formed under general conditions. Slow switching appears as a result 
of the system's approach to the saddle connections. It is also argued that such saddle connections 
easy to arise near the bifurcation point where the state of perfect synchrony loses stability. We 
develop an asymptotic theory to reduce the model into a simpler form, with which an analytical 
study of cluster states becomes possible. 
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C3 . I. INTRODUCTION 

CZ5 . 

Studies on collective motion of coupled oscillators have attracted considerable attention over the last few 
decades P, 0, 0- It is commonly seen that a population of autonomous elements performs certain biological func- 
tions by behaving collectively!^. It has in fact been pointed out that collective motion is crucial to information 
processing and transmission in living organisms 0. 

In the brain, the neurons are exclusively coupled through chemical synapses, i.e., the neurons communicate by pulses 
of transmitter |(|. Chemical synapses commonly form dense and complex networks. For mathematical modeling of 
neuronal networks, homogeneous all-to-all (or, global) coupling is often adopted. Although the global coupling may 
be a little too idealistic, the corresponding networks share a lot of properties in common with systems with complex 
and dense networks. 

In the present paper, we consider a population of neural oscillators with delayed, all-to-all pulse-coupling. The 
oscillator we use is called the leaky integrate-and-fire (LIF) model. There are a large amount of papers concerning 
LIF in physics and neuroscience, e.g., see 0,11, This is because LIF is a quite simple model still capturing some 
1 essential characteristics of neuronal dynamics, i.e., it represents an integrator with relaxation, and resets after it fires. 
Though our population model is commonly used, (e.g., see [Tol|'). its collective dynamics does not seem to have been 
l/S ■ studied so carefully. We are particularly concerned with peculiar collective dynamics called slow switching [Til IT^|. 

The study of collective dynamics in the original form of the model is not easy to handle because the coupling involves 
a long term memory. We thus develop an asymptotic theory and reduce our model into a form without memory, by 
' which an analytical study of collective dynamics becomes possible. 

Eh ; II. MODEL 

The population model we consider consists of N identical elements with delayed, all-to-all pulse- coupling. The 
dynamics of each elements is described by a single variable (i = 1, 2, . . . , N) which corresponds to the membrane 
J"^ | potential of a neuron. The equation for Vi is given by 



±v i (t)=a-v i + ^(b-v i )E(t). (1) 

The parameter a is the so-called resting potential to which Vi relaxes when the coupling is absent. It is assumed 
that when u,; reaches a threshold value which is set to 1, it is instantaneously reset to zero. This event is interpreted 
as a spiking event. The dynamics is thus called LIF. When a neuron spikes, it emits a pulse toward each neuron 
coupled to it, and the latter receives the pulse with some delay called a synaptic delay. The coupling is assumed to 
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be homogeneous and all-to-all, so that its effect can be represented by one global variable E, given by 

N 

£(*) = E E e(t-tf ies -r). (2) 

J=l spikes 

Here, ^ plkes represents a series of times at which the j-th neuron spikes and X^spikes denotes a summation over the 
series of such spikes; r is the synaptic delay, and e(t) is a pulse function, given by 

e(t) = -^(e- at -e-^)Q(t), (3) 
p — a 

where Q(t) is the Heaviside function: a and /? are constants satisfying (3 > a. In the limit /3 — > a, e(i) becomes a 2 te~ at , 
which is called the alpha function 6J. 6 is called the reversal potential to which Vi relaxes when E(t) is positive, i.e., 
while the neuron receives the pulses. K is a positive constant characterizing the strength of the coupling. The 
coupling assumed above is characteristic to the synaptic coupling. The coupling become excitatory (EPSP) if b > 1, 
and inhibitory (IPSP) if b < 0. 

If a < 1, LIF becomes an excitable neuron, while if a > 1, it repeats periodic spikes, namely, it represents a neural 
oscillator. We assume a > 1 throughout the present paper, and call each element an oscillator. Then, we can define 
a variable ipi varying smoothly with time, which turns out useful in the following discussion. We call ipi the phase of 
the i-th oscillator, and define it by 

f— =i*(— y (4) 

Jo a-v \a~ViJ 
which varies between and the intrinsic period of oscillation T given by 

T-uUrY (5) 



Note that ipi satisfies dipi/dt = 1 in the absence of coupling. 



III. NUMERICAL RESULTS 



By numerically integrating our model under random initial distributions of Vi, we find various types of collective 
behavior. Among them, we are particularly interested in the slow switching phenomenon, which can arise when b > a 
and N > 4. As displayed in Fig. ^ the whole population, which was initially distributed almost uniformly, splits into 
two subpopulations, each of which converges almost to a point cluster. However, after some time the phase-advanced 
cluster starts to scatter. Then, this scattered group starts to converge again as it comes behind the preexisting cluster. 
In this way, the preexisting cluster becomes a phase-advanced cluster. After some time, again, this phase-advanced 
cluster begins to scatter, and a similar process repeats again and again. In other words, the system switches back 
and forth between a pair of two-cluster states. For larger times, the system comes closer to each of well-defined 
two-cluster states and stays near the state longer. Theoretically, these switchings repeats indefinitely, although in 
numerical integrations the system converges at one of the two-cluster states in a finite time and stops switching due 
to numerical round-off errors p^. 

The slow switching phenomenon occurs within a broad range of parameter values provided that K is small, and the 
time constants a -1 , /3 _1 and r are small compared with T. For larger a^ 1 , /3 _1 and r, the slow switching phenomenon 
becomes less probable, and the appearance of steady multi-cluster states becomes more probable instead. For 6 < a, 
we find no two-cluster states involving slow switching, while steady multi-cluster states are observed in most cases. 
The corresponding phase diagram will be presented in Sec. IVIII (see Fig. Q). 



IV. WEAK COUPLING LIMIT 



Our model given by Eq. is relatively simple, still it would be difficult to get some insight into its collective 
dynamics analytically. Fortunately, however, our main results given in Sec. IIIII do not change qualitatively in the 
weak coupling limit, i.e., K — » 0. In this limit, our model is reduced to a much simpler form with which we can study 
the existence and stability of various cluster states analytically. Derivation of the reduced model is given as follows. 
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Substituting m = a(l — e into Eq. (JIJ, we obtain 



where 



f t Mt) = i + £ E E ^M* - *f ikes - ( 6 ) 

j — 1 spikes 



Z(V><) = — e*< + 1. (7) 



It is convenient in the following calculation to redefine Z as a T-periodic function, or, + nT) — Z(%pi) [n — 

±1, ±2, • • ■ ). Note that sudden drop of Z(x) at x — is due to our rule employed, i.e., the membrane potential is 
instantaneously reset at Vi = 1. We also define a residual phase \&i by 

= - *• (8) 

Substituting Eq. (jHJ into Eq. (JHJ, we obtain 

H K N 

J t ^) = Jjz2 E m+t)e(t-tf^ s -r). (9) 

J — 1 spikes 

We now assume that K is sufficiently small so that the r.h.s of Eq. © is sufficiently smaller than the intrinsic frequency 
T -1 . This allow us to make averaging of the r.h.s of Eq. @ over the period T. The zeroth order approximation with 
respect to the smallness of K, which corresponds to the free oscillations, is given by 

*i(t) = const. (10) 

and 

if kos = t,-nT, (n = 0,1,2,...), (11) 

where tj is the latest time at which the j'-th neuron spikes. In the first order approximation, we may time-average 
Eq. © over the range between t — T and t using Eqs. (|10|l and (fill) : 

d , , , KAl '* °° 



dt 9i ® = Ifil^W f Il Z (Mt) + X)e(\-t 3 +nT-T)dX (12) 

j=l • !t - T n=0 

= jvEy/ + + r + ^')e(A')dA' (13) 



f + d4) 



3=1 



where 



r( x ) = ^L{ ff _ ^, T ^)}, (15) 

p — a 



1 r 00 

Hu,t{x) = — exp[(x + A) mod T] exp[— a\]d\ 

2 Jo 

(e T - 1) exp[a(x mod T)] - (e aT - 1) exp[a; mod T] 



T(l-a)(e QT -l) 



(16) 



Note that r(x) and H a ^{x) are T-periodic functions. Figure [21 illustrates a typical shape of the coupling function 
given by Eq. 115|) . Furthermore, using the identity 

*j(*i)=^(*i)-*i = -*i. (17) 
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and the zeroth order approximation ^j(tj) — ^j(t), we may replace tj by — ^j(t) in Eq. (|14|l in the first order 
approximation. Thus, we finally obtain 

d K' N 

-Ut) = u + £ r (V^) - + t)> ( 18 ) 

3=1 

where w = 1 + i^/T and K' = K (b — a) /a. Equation l|18(l is the standard form of the phase model. Note that the error 
involved in Eq. may look to diverge as n — > oo, still the final error vanishes in the first order approximation due 
to the decay of e(t). It should be noted that the reduced model is free from memory effects, but the effect of delay 
has been converted to a phase shift in the coupling function. Similar form of the phase model is generally obtained 
in delayed coupled oscillators when the coupling is sufficiently weak^^- Hereafter, we ignore the degree of freedom 
associated with the dynamics of the center of mass (or, mean phase) which can be decoupled in the phase model. 

Important parameters of our phase model given by Eq. (|18|l with Eq. I|15|) are T,a,(3,T and the sign of K' (i.e., 
the sign of b — a). The reason is the following. We may take \K'\ = 1 by properly rescaling of t and u, while its 
sign is crucial because the local stability of any solution depends on it. u> gives the frequency of steady rotation of 
the whole system, which is irrelevant to collective dynamics. We choose T as an independent parameter by which a 
becomes dependent through Eq. JSJ. It is remarkable that our coupling function is independent of b. In fact, change 
in b causes no qualitative change in our result as far as the sign of b — a remains the same. Interestingly, even if we 
replace the term b — Vi by a constant c in Eq. (Q, i.e., 

d Kc 

- Vi (t) = a-Vi + —E(t), (19) 

then we can reduce this model similarly and obtain the same coupling function as in Ea. (|15|l . We have checked 
that Eq. (|19|) actually reproduces qualitatively the same results as those given in Sec. IIIII In that case, negative c 
corresponds to the case b < a in Eq. (JIJ. 

In the following section, we assume b > a and (3 — > a unless stated otherwise. 



V. TWO-OSCILLATOR SYSTEM 



In this section, we study a two-oscillator system, or, N — 2. Although the two-oscillator system is not directly 
related to the main subject of the present paper, one may learn some basic properties of our phase model from this 
simple case. Defining A = ?/>i — "02 > we obtain 

^ = ^(r(A + r) - T(-A + r)) = G T (A). (20) 

Phase locking solutions are obtained by putting dA/dt = 0, and the associated eigenvalues are given by dG T /dA. 
Figure |31 shows a bifurcation diagram of the phase locking solutions, in which we take r as a control parameter. We 
find that for small r the trivial solutions A = (in-phase locking) and T/2 (anti-phase locking) are unstable, while 
there are a pair of stable branches of non-trivial solutions. The point r = is close to the bifurcation point where 
the in-phase state loses stability. The bifurcation occurs at r = t c , where r c corresponds to the minimum of T(x) 
(see Fig- EJ) - Because r c is negative, the in-phase state cannot be stable for small or vanishing delays (while it can be 
stable for delays comparable to T due to the T-periodic nature of our phase model). r c is extremely small, which is 
due to the sudden drop of Z(x) at x = and the particular rule employed in our model, i.e., a neuron is assumed to 
spike and reset simultaneously. The width of the stable branches of the trivial solutions is the same as that of the 
decreasing part of T(x). Owing to the peculiar shape of Z(x), the width is of the same order as the width of e(t), 
which is 0(a^ 1 ). The stability of the in-phase state is identical with that of the state of perfect synchrony. 

In terms of the original model, we now present a qualitative interpretation of why the in-phase locking state is 
unstable for small or vanishing delays. We consider the dynamics of two neurons which are initially very close in 
phase. The effect of a pulse on the phase ipt is larger for smaller dvi/dt. dvi/dt is monotonously decreasing except 
when it is reset (which reflects on the property of Z(x) that it is increasing except x = 0). Thus, the neuron with 
larger Vi makes a larger jump in phase when it receives a pulse, by which the phase difference between the two neurons 
becomes larger when they receive a pulse. On the other hand, the situation becomes different if two neurons lie before 
and after the resetting point, i.e., if the phase-advanced neuron has smaller Wj. In that case, the phase difference 
becomes smaller when they receive a pulse. According to our dynamical rule, however, resetting and spiking occur 
simultaneously, so that they receive pulses when the phase-advanced neuron has larger Uj. Therefore, the in-phase 
state becomes inevitably unstable even without delay. If we want to obtain a stable in-phase state for small delays, 
we should employ a rule such that a neuron spikes before it is reset, which would be more physiologically plausible 
than the rule employed here. 
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VI. LOCAL STABILITY ANALYSIS FOR A LARGE POPULATION 

The trivial in-phase solution and the non-trivial solutions of the two-oscillator system correspond respectively to 
the state of perfect synchrony and two-cluster states when we go over to a large population. In this section, we study 
local stability of the two-cluster states. Although non-trivial solutions are stable for small or vanishing delays in the 
two-oscillator system, the corresponding two-cluster states turn out unstable. 

We consider a steadily oscillating two-cluster state in which the two clusters consist of Np and N{\ — p) oscillators, 
respectively. The oscillators inside each cluster are completely phase-synchronized, and the phase-difference between 
the clusters is constant in time, which is denoted by A. From Eq. I|18|l . the phase difference obeys the equation 

^=K' {(2p - l)r(r) + (1 - p)T(A + r) + P r(-A + r)} . (21) 

When A is constant, we obtain a relation between p and A as 

(2p - l)r(r) + (1 - p)T(A + t)+ pT(-A + r) = 0. (22) 

We designate a two-cluster state satisfying Eq. lt^|l as (p, A). The eigenvalues of the stability matrix are calculated 

as 

X 1 = K'{pT'(r) + (l-p)r'(A + r)}, (23) 



A 2 = K'{(1 - p)T'(r) + pT'(-A + r)}, (24) 



A 3 = K'{{1 - p)T'{A + r) + pT'{-A + r)}, (25) 

where r'(a;) means (d/dx)T(x). The multiplicities of Ai, A2 and A3 are Np— l,iV(l — p) — 1 and 1, respectively. Ai and 
A2 correspond to fluctuations in phase of the two oscillators inside the phase- advanced and phase-retarded cluster, 
respectively. A3 corresponds to fluctuations in the phase difference A between the clusters. 

Substituting Eq. (JT3J) into Eq. O, we obtain a relation between p and A, the corresponding curve being shown 
in Fig. Ola). By using this relation, the eigenvalues of (p, A) can be obtained, which are displayed in Fig. 0Jb) as 
a function of A. It is found that no two-cluster states are stable, and there is a set of (p, A) for which Ai > and 
A2, A 3 < 0. Positive Ai means that the two-cluster state is unstable with respect to perturbations inside a phase- 
advanced cluster. On the other hand, A2, A3 < means that it is stable against perturbations inside a phase-retarded 
cluster as far as the perfect phase-synchrony of the phase-advanced cluster is maintained. Within a certain range of 
p, there are pairs of two-cluster states represented by (p, A) and (p, A') with the same stability property, and they 
appear as the solid lines in Fig^a). In the next section, we explain how a heteroclinic loop between the (p, A) and 
(p, A') is stably formed in our model. 

Similarly to the discussion in Sec.0 the stability property mentioned above can also be understood in terms of the 
original model. Every neuron inside the phase-advanced cluster always receives pulses when its membrane potential 
is increasing. Then, the phase-difference between two neurons inside the cluster, if any, always increases, so that the 
phase- advanced cluster is inevitably unstable. On the other hand, the neurons inside the phase-retarded cluster can 
receive pulses (emitted by the phase- advanced cluster) during their resetting. Then, the phase-differences between 
neurons inside the phase-retarded cluster, if any, become smaller, so that the phase-retarded cluster can be stable. 

VII. HETEROCLINIC LOOP 

We first note that there is a particular symmetry of our model which turns out crucial to the persistent formation 
of the heteroclinic loop. The symmetry is given by 

= for any i and j. (26) 

»<(<)=«<(*) 

Due to this symmetry, the units which have the same membrane potential at a given time behave identically thereafter. 
In other words, once a point cluster is formed, it remain a point cluster forever. 
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We assume that a pair of two-cluster states (called A and B) exists and has the same stability property as that 
discussed in Sec. I VII i.e., the phase-advanced cluster is unstable, and the phase-retarded clusters is stable. Suppose 
that our system is in a two-cluster state A initially. When the oscillators inside the phase-advanced cluster are 
perturbed while the phase-retarded cluster is kept unperturbed [see Fig.|Sfa)], the former begins to disintegrate while 
the latter remains a point cluster. Then, the group of dispersed oscillators and the point cluster coexist in the system 
[see Fig. Elb)]. We know, however, that in the presence of a point cluster, there exists a stable two-cluster state in 
which the existing point cluster is phase-advanced. From this fact, the dispersed oscillators are expected to converge 
to form a point cluster coming behind the preexisting point cluster. In this way, the system relaxes to another two- 
cluster state B [see Fig.|Jc)]. From the above statement, it is implied that in our high-dimensional phase space there 
exists a saddle connection from the state A to the state B. The existence of a saddle connection from the state B to 
the state A can be argued similarly. A heteroclinic loop is thus formed between the pair of the two cluster states A 
and B. 

In terms of the phase model, the above argument can be reinterpreted in a little more precise language. In the 
phase model given by Eq. . a symmetry property similar to Eq. I|26|) also holds: 



d 
dt 



= for any i and j. (27) 

V>i(t)=<M*) 



Our argument will be based on the following assumptions: 

(a) (p, A) with Ai > and A2, A3 < exits, 

(b) (p, A') with A' 2 > and X[, A 3 < exits, 

where we define A > and A' < 0, and Ai and A- (i = 1, 2, 3) are the eigenvalues of (p, A) and (p, A'), respectively. 
Note that if p — 0.5, the two clusters in question are identical, or A' = A, so that (a) and (b) are identical. Figure 
[5] illustrates a schematic presentation of the N — 1 dimensional phase space structure, where we ignore the degree of 
freedom associated with the dynamics of the center of mass. E a and E r are identical with the subspaces where the 
phase-advance and phase-retarded clusters of (p, A) remain point clusters, respectively. By considering the direction 
of eigenvectors, one can easily confirm that E a and E r are identical with the stable subspaces of (p, A) and (p, A'), 
respectively. Furthermore, because E a and E r are invariant subspaces due to the symmetry given by Eq. J57J), 
(p, A) and (p, A') are attractors within E a and E r , respectively. Thus, a heteroclinic loop between (p, A) and (p, A') 
should necessarily exits. The saddle connections in question are stably formed through the invariant subspaces which 
exist for the symmetry of equations of motion given by Eq. I|27|) . The heteroclinic loop is thus robust against small 
perturbations to the system unless the symmetry is broken. 

Whether the resulting heteroclinic loop is attracting or not depends on the following quantity: 

It was argued in Ref. 01 that if 7 > 1, the system can approach the heteroclinic loop and come to move along it. 
In that case, the time interval during which the system is trapped to in the vicinity of one of the two-cluster states 
increases exponentially with time. Substituting the eigenvalues obtained from Eqs. (|23|) and lt2^|) using Eq. (|T3|) into 
Eq. I|28f) . we find that the heteroclinic loops within a certain range of p are in fact attracting for small a , /3 _1 and 
t. Phase diagrams of the heteroclinic loops and symmetric multi-cluster states is shown in Fig.[7| where we choose r 
as a control parameter (see appendix for the stability of the symmetric multi-cluster states) . 



VIII. NEAR THE BIFURCATION POINT 



In this section, we concentrate on the vicinity of the bifurcation point where the state of perfect synchrony loses 
stability. As noted in Sec. the bifurcation occurs at r = r c . Then, for small x — r c , the coupling function can be 
expanded as 

T(x) = co + c 2 (x - r c ) 2 - c 3 (x - r c f + 0((x - r c ) 4 ). (29) 

Suppose that c 2 and C3 are positive. We further put C3 = 1 by properly rescaling K' in Eq. Q18[l. In order to find 
possible two-cluster states, we solve Eq. lEffi using Eq. (|2"§|) . We then obtain three solutions for A as a function of p 
and r. One is the trivial solution A = (the perfect synchrony), and the others are given by 



A = (l-2p)(<»-3f) ± /(l-2^(c 2 -3f)3 + ^ 
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where f = r — r c . Note that the expression above using the approximate T given by Eq. (|29(l is valid only for small A, 
which is actually the case if p is close to 1/2 and f is small. Substituting the expressions in Eq. (|3*U)l into Eqs. (|23l - (|2l))l . 
we obtain eigenvalues associated with the two-cluster states. The resulting bifurcation diagram for given p is shown 
in Fig. [S] The solid and broken lines give the branches of negative and positive A3 , respectively. Two solid branches 
exist for r > 0, which are represented by (p, A) and (p, A') with A > and A' < 0. One can easily confirm that the 
eigenvalues of these states satisfy Ai, X' 2 > and A2, A3, A' 1; A 3 < for arbitrary p and small f , which agree with the 
condition for the existence of a heteroclinic loop. The quantity 7 defined by Eq. (|28[) can also be calculated and turns 
out to be larger than 1. Thus, all the local stability conditions for the existence of an attracting heteroclinic loop are 
generally satisfied just above the bifurcation point provided C3 > 0. 

It is also possible that a heteroclinic loop is formed when c 3 < 0. In that case, it is expected to arise subcritically, 
so that both the heteroclinic loop and the state of perfect synchrony may be stable over some region of negative f . In 
fact, we found that such bistability arises in a population of the Morris-Leccar oscillators with the same coupling 
form as in Eq. JQ), and an analysis by means of the phase dynamics actually shows that C3 is negative. To confirm 
the corresponding bifurcation structure, we have to consider higher orders of x — t c in the coupling function. The 
details of this issue are omitted here. 



IX. CONCLUSIONS AND DISCUSSION 



We have discussed the slow switching phenomenon in a population of delayed pulse-coupled oscillators. We found 
that the phenomenon is caused by the formation of an attracting heteroclinic loop between a pair of two-cluster states. 
A particular stability property of the two-cluster states and a certain symmetry of our model are responsible for its 
formation. Our original model given by Eq. is reduced to the standard phase model in the weak coupling limit, 
by which we succeeded in studying the stability of the two-cluster states analytically, and confirming the structure of 
the heteroclinic loop. It was also argued that under the mild condition of the coupling function all the local stability 
conditions for the existence of an attracting heteroclinic loop are generally satisfied just above the bifurcation point. 

The physical mechanism of the formation of a heteroclinic loop we describe in Sec. IVHI does not depend on the 
nature of elements (e.g., phase oscillator, limit-cycle oscillator, excitable elements, chaotic elements) and couplings 
(e.g., diffusive coupling, pulse coupling). It is expected, therefore, that a heteroclinic loop arises in a wide class of 
models of coupled elements. 
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APPENDIX A 

According to Ref.|l3l|. we summarize here the existence and the stability analysis of symmetric multi-cluster states 
in the phase model given by Eq. I|18(l . In the symmetric n-cluster state, it is assumed that each cluster consists of N/n 
oscillators. We denote the phase of cluster k as <pk {k = 0, 1, . . . , n — 1). There always exists the following solution: 

Tk 

<f>k=ttt+ — , (Al) 
n 



K' (Tk 



with 

^ = -£ r (-+-h (A2) 

which corresponds to the state in which the n clusters are equally separated in phase and rotate at a constant frequency 
il. The associated eigenvalues are calculated as 

A intra =^yW— +T ) ) (A3) 



8 



| number density 
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FIG. 1: Slow switching phenomenon viewed through the number density of the oscillators as a function of phase. In order to 
get a better view, we work with a comoving frame of reference. The parameter values are a — 1.03 (T ~ 3.5), b = 2.0, aT 1 — 
0.3, /3 -> q,t = 0.2, K = 0.1 and N = 100. 



K' ^4 , /Tfc 



A fnt cr = — E r ' hr + T M 1 - exp[-*Tfcp/n]). (A4) 
11 k=o v " 



Aintra is a intra-cluster eigenvalue with multiplicity of N — n. Af nter (p=l,. . . ,n-l) are associated with inter-cluster 
fluctuations. If all of these eigenvalues have negative real part, the symmetric n-cluster state is stable. 
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FIG. 2: The solid line shows the coupling function T(x) for a = 1.05 (T ~ 3.0), a 1 = 0.2 and f3 — > a. The minimum appears 
at cc = t c which is a small negative. For comparison, the shape of Z(x) is also displayed with the broken line. 
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FIG. 3: Bifurcation diagram of a two-oscillator system. Solid and dotted lines respectively represent stable and unstable 
branches, where b > a is assumed. The stability property is reversed if b < a. 
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FIG. 4: (a)Relation between p and A associated with two-cluster states, (b) Eigenvalues of two-cluster states as a function of 
A. In (a), the solid and dotted lines correspond to the two-cluster state of negative and positive A3, respectively. 



10 




FIG. 5: Schematic representation of a saddle connection between a pair of two cluster states, starting with the two-cluster state 
A (a) ending up with the other two-cluster state B (c). 




FIG. 6: Schematic representation of the structure of a heteroclinic loop, (p, A) and (p, A') become attractors in the invariant 
subspaces E a and E r , respectively. 
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FIG. 7: Phase diagrams of cluster states, where r is chosen as a control parameter. The parameter values are the same as 
in Fig. with (a)& > a and (b)6 < a, respectively. For given p and r inside the gray region, 7 is larger than one, i.e., the 
heteroclinic loop between (p, A) and (p, —A') is attracting. Each rectangle placed at n indicates the region of r within which 
the symmetric n-cluster state is stable. In (b), stable symmetric n-cluster states with n > 10 also exist (not shown). 
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FIG. 8: Bifurcation diagram around r = r c . A heteroclinic loop is formed between a pair of the solid branches for r > r c . 



